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Abstract. We benchmark many-body perturbation theory by studying neutral, 
as well as non-neutral, excitations of finite lattice systems. The neutral excitation 
spectra are obtained by time-propagating the Kadanoff-Baym equations in the 
Hartree-Fock and second Born approximations. Our method is equivalent to 
solving the Bethe-Salpeter equation with a high-level kernel while respecting self- 
consistently, which guarantees the fulfillment of a frequency sum rule. As a result, 
we find that a time-local method, such as Hartree-Fock, can give incomplete 
spectra, while already the second Born, which is the simplest time-nonlocal 
approximation, reproduces well most of the additional excitations, which are 
characterized as double-excitations. 



PACS numbers: 31.10.+z,71.10.-w, 31.15.xm 

1. Introduction 

We recently developed a non-equilibrium many-body method, based on time- 
propagation of the Kadanoff-Baym equations (KBE), to study transient dynamics in 
quantum transport, and found that electronic correlations beyond mean-field strongly 
affect the relevant non-equilibrium properties [TJ [2]. The method is approximate, 
because only some selected classes of perturbative terms are accounted for, and 
therefore an estimate of the quality of the results would be desirable. However, 
as benchmark results are scarcely available for such open systems, we need to 
devise simpler experiments to test our method. Some finite, exactly solvable lattice 
systems, which typically model molecular devices and whose properties thus determine 
the transport characteristics, form an appropriate test platform. A many-body 
approximation which is unable to reproduce the properties of such an isolated system 
correctly, has no prerequisites to fare well for a molecular device between conducting 
leads. 

The Kadanoff-Baym equations are equations of motion for the one-body Green's 
function G, where many-body effects are incorporated into the so-called self-energy 
E. The latter is a functional of the Green's function, and the main task of many-body 
perturbation theory (MBPT) is to find good approximations for this quantity. We are, 
in particular, interested in the so-called conserving approximations for the self-energy, 
which are vital in quantum transport [3l |4] as they guarantee the fulfillment of the 
conservation laws for the particle number, total momentum and energy [5l|6]. Some 
of these approximations have already been tested in finite lattice systems [U |8] . These 
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reports have highlighted strongly perturbed systems, in which serious issues, such as 
artificial damping of the dynamics, were observed and deemed correlation-induced as 
they were not present in a mean-field approximation. The purpose of this paper is to 
extend these investigations to the linear response regime. 

Within this regime, the key quantity is the response function 6G/Sv which 
describes the change in the Green's function due to a weak external field v, and 
gives a direct access to the neutral excitation spectra. This quantity satisfies an 
integral equation known as the Bethe-Salpeter equation (BSE) [9], whose integral 
kernel is given by the functional derivative 6T,/6G. Using this equation to calculate the 
excitation spectra has proven to be computationally challenging and, in practice, often 
requires a number of additional approximations, such as neglect of self-consistency, 
kernel diagrams, and/or frequency-dependence. Instead, we obtain the response 
function by time-propagation of the Green's function which does not require any of the 
above mentioned approximations [inilll]- Moreover, such obtained excitation spectra 
automatically satisfy a frequency sum rule known as the f-sum rule, which acts as an 
important consistency relation for the response function. 

The spectral properties of neutral, as well as non-neutral, excitations depend 
strongly on the temporal structure of the self-energy which can be either time-local 
or -nonlocal. Whereas a time-local scheme leads to merely renormalized one-particle 
excitations, a time-nonlocal one can reproduce more complicated spectral structure. 
Such excitations involving double- or many-particle character are, in fact, dominant for 
some potential materials for realistic devices, such as organic polymers (polyenes) |12) . 
making time-nonlocal approximations an interesting research topic. This represents 
a major challenge for the Bethe-Salpeter approach as time-nonlocality translates to 
an integral kernel which depends on three frequencies. Although some state-of-the- 
art calculations have been recently conducted in finite systems with a time-nonlocal 
kernel [13l [M], these results have become at the expense of self-consistency which 
is a requirement for fulfillment of conservation laws, and therefore indispensable in 
a time-dependent transport process. Moreover, a non-self-consistent approximation 
can lead to a gross violation of the f-sum rule [TOl [H]. On the other hand, our 
calculations, conducted within the time-local Hartree-Fock (HF) and time-nonlocal 
second Born (2nd Born, 2B) approximations [HI [16], do not sacrifice self-consistency, 
remain computationally tractable, and yield excitation spectra satisfying the f-sum 
rule. 

The paper is organized as follows. We start with an introductory section [2] in 
which we summarize the non-equilibrium many-body formalism, describe our means 
to approach the neutral excitation spectra, and introduce, as well as analyze, our 
many-body approximations. In section |3l we introduce the test environment which 
comprises two simple lattice systems. Our numerical results, in particular, on addition, 
removal, and excitation spectra, as well as some related discussion are contained 
in section U] We summarize our work and present some conclusions in section [5] 
Finally, [Appendix A contains an extension of the proof the f-sum rule to lattice 
systems. 
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2. Response functions in non-equilibrium many-body theory 

2.1. Non- equilibrium Kadanoff-Baym approach 

The Green's functions are reduced quantities which act as probes to physical 
observables, while containing only a minimal amount of excess information. The 
non-equilibrium one-body Green's function is defined as 

Gij{z,z') = -i{%,CHAz)cH,j{z')) , (1) 
where (...) denotes an ensemble average with respect to the grand canonical density 
operator, and ^^■'^(z) a Heisenberg picture operator which annihilates (creates) a 
particle from (to) a single-particle quantum state i. The contour-ordering operator Ty 
arranges the contour times z and z' along the extended Keldysh contour 7 [TT^ shown 
in figure [U for further details see [HI . 

The contour-ordered Green's function, although a powerful formal device, does 
not represent directly any observable. It can, however, be reduced to a set of more 
intuitive and physical objects by exposing all possible time-orderings. The one-body 
Green's function, being a two-time object, can then be written as 

G{z, z') ^ 9{z, z')G> {z, z') + e{z' , z)G<{z, z') , (2) 

where the boldfaced symbols denote matrices in the basis spanned by the single- 
particle quantum states. The greater and lesser Keldysh components, which are 
defined as 

G>{z,z') = -i{cH4z)c'H^jiz')) , (3a) 

G<(z,z') = i{c'^ .j{z')cH4z)) , (3&) 

describe the motion of a particle and a hole, respectively, in the many-particle system. 
The core of the non-equilibrium theory lies in its equations of motion, the first of 
which can be written as 

{id, - h{z))Giz, z') = S{z, z') + J dz S[G](2, z)G{z, z') , (4) 

7 

while the second, adjoint equation can be obtained by differentiation with respect 
to the second time-argument. The symbol h{z) denotes the one-body part of 
the Hamiltonian which can, in general, be time-dependent. Using the Langreth 



Figure 1: The extended Keldysh contour starting from to and ending to to — mif3 
consists of an equilibrium, imaginary track and a two non-equilibrium, real-time 
branches. The zero-temperature limit /? — > 00, where /3 is the inverse temperature, is 
implied when finite systems are considered. 
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rules [IH], the equations of motion can be split into a set of equations for the 
different Keldysh components referred as the Kadanoff-Baym equations [HI HH] ■ The 
self-energy S is an effective, nonlocal potential, which accounts for all the many-body 
effects, and allows a systematic, beyond order-by-order perturbation expansion. A 
Green's function obtained from an approximate self-energy 

where is the Baym- functional [6], obeys the conservation laws for the particle 

number, total momentum and energy provided that it satisfies the equation of 
motion or in other words is solved self-consistently fSl, 

The nonlocal structure of the one-body Green's function ensures that all one- 
body observables which are accessible in the equilibrium/zero-temperature theory can 
be now calculated in non-equilibrium settings from the knowledge of the greater and 
lesser components. In particular, all time-local observables are related to the one-body 
reduced density matrix (IRDM) which is given by 

7(i)^-iG<(t,t), (6) 

where t is a real-time argument located on the horizontal track of the Keldysh contour. 
On the other hand, time-nonlocality allows access to information on particle number 
changing processes contained in the spectral function, which is defined as 

A{t,t')^i[G>it,t')~G<it,t')] , (7) 

and whose pole structure in the frequency-domain gives the particle addition and 
removal energies. Another virtue of time-nonlocality, is the possibility to access some 
time-local two-body observables, including the total energy which can be written in 
terms of the Galitski-Migdal (GM) functional as 

E{t)^-'-tY\(dt + h)G<{t,t')\ . (8) 



2 



t'=t 



where dt is proportional to the identity matrix, and tr A = J2i is the regular 
algebraic trace. However, the main advantage of the non-equilibrium theory addressed 
here, is the additional possibility to calculate some time-nonlocal two-body observables 
from the sole knowledge of the one-body Green's function. We focus on the real-time, 
retarded response function, which is defined as the functional derivative 



5vik{t') 



(9) 

with respect to an external perturbation Vij{t), and has the general structure of a 
retarded function given by 

xlu it, t') = e{t-t') (x> (t, t') - x< (t, t')), (10) 

where the greater and lesser components are defined, respectively, in terms of the 
one-body reduced density matrix operator 'jH.ij 

(t) = c'jj .{t)cH4t), as 
X>,kiit,t') =-i(7ff,,(07H,M(i')), (11a) 

xf,,kii^,z') ^ -i{iHMit')iHMt)) ■ im 

The retarded response function, therefore, describes the motion of a particle-hole 
pair, or in other words an excitation, in the many-particle system. This quantity is 
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connected to the one-body Green's function through the hnear response relation which 
can be written as 



oo 




where S^ij (t) is the first-order variation of the one-body reduced density matrix. The 
lowest order response of the system to an external perturbation is, therefore, uniquely 
determined by the real-time, retarded response function, whose pole structure in 
the frequency-domain gives the neutral excitation energies. In the non-equilibrium 
approach, both the IRDM and perturbation are known, and therefore the response 
function is readily obtained by a proper choice of perturbation, and subsequent 
inversion of equation (IT^ pUl ITT] . 



2.2. Generalized response function 

The quality of the one-body Green's function and, because of equation also of 

the retarded response function is clearly determined by the self-energy approximation. 
Here we write down an explicit relation between the self-energy and retarded response 
function. We start by introducing the generalized, contour-ordered response function, 
which is defined as 



(13) 

v=0 



6vik{z") 

and in terms of which, the greater and lesser components of the real-time response 
function, defined in equation p.ip . are given by 

(14) 

-t± ,z'—t' 



Xfj,kli*,t') = -iLjj-fc;(zz + ,z') 



where z+ = z + ri with the limit 77 — 0+ implied after contour-ordering. The contour- 
times z — t± are located on the lower/upper horizontal branch of the Keldysh contour 
shown in figure [TJ and correspond to a real-time t. The generalized response function 
can be straightforwardly shown 9 to obey the so-called Bethe-Salpeter equation which 
is given by 

Lij^ki{zz' , z") = Gti{z, z")Gkjiz" , z') 

+ d{ziZ2Z3Z4) Gip{z, Zi)Gqj{z2, z')Kpq^sr{ziZ2, Z4Z3)Lrs,kl(^3^4:, z") , (15) 



pqrs ^ 

where we defined the four-point Bethe-Salpeter kernel as the functional derivative 

Figure [2] contains a diagrammatic representation of this equation, and illustrates the 
role of the kernel as a source of an infinite perturbation series. We can conclude that 
any retarded response function obtained by inversion of equation (IT2|) . and relating 
to a one-body Green's function for an approximate self-energy E, corresponds to a 
solution of the Bethe-Salpeter equation with the kernel 5Y,/5G. The correspondence 
further implies that a non-equilibrium approach already with a relative simple self- 
energy, accounts to a quite sophisticated approximation for the response functions. 
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Figure 2: The Bethe-Salpeter equation for the generalized response function whose 
particular time-ordering Lij^ki{zz~^ , z') known as the particle-hole propagator yields 
the retarded response function. The dashed lines represent possible connections, while 
double lines denote dressed Green's functions. 



2.3. Approximate Bethe-Salpeter kernels 

The time-propagation of a Green's function with an approximate self-energy S 
automatically leads to a retarded response function related to an approximate kernel 
6Yi/5G which, for example, guarantees the obedience of the f-sum rule, as proven 
in [Appendix A[ Such consistency ensures that each self-energy diagram containing n 
Green's function lines, gives rise to n inequivalent kernel diagrams which contribute to 
the response function. Our two approximate kernels are given by functional derivatives 
of the Hartree-Fock and 2nd Born self-energies, and represent a time-local mean-field 
and a nonlocal correlated approximation, respectively. 

We can view our generalized response functions diagrammatically by inserting 
the approximate kernels into the Bethe-Salpeter equation shown in figure [2] Guided 
by this recipe, we easily deduce that the Hartree-Fock kernel shown in figure [3a| yields 
a response function given by a simple bubbles-and-ladders series. Whereas the second 
Born kernel, which is given in figure [3bl results into a considerably more complicated 
response function. Note that, our approximations, whether related to a self-energy 
or a kernel, are known by the name of the self-energy without any extensions. For 



Hartrco (H) Fock (F) 



HF = > ' " — ^ l^HF 



Khf = 



(a) Hartree-Fock: self-energy Shf and Bethe-Salpeter kernel Khf 
S'^hf/Ghf- 



bubble (B) exchange (X) 




Bi— B3 Xi— X3 

(b) 2nd Born: self-energy J^2B and Betlie-Salpcter kernel K2B = SY^2B jG^B- 



Figure 3: Self-energy diagrams give rise to kernel diagrams, for example bubble (B) 
gives diagrams Bi_3. Wiggly lines denote the Coulomb interaction, and double lines 
dressed Green's functions. 
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example, our Hartree-Fock approximation for the kernel shown in figure |3a] is often 
explicitly referred to as the time-dependent Hartree-Fock approximations. 

We can also assign a meaning to each diagram in these perturbation expansions by 
considering a particular time-ordering which, in our case, is the ordering given by the 
particle-hole propagator Lij^ki{zz~^ , z'). We interpret our kernel diagrams as follows. 
The first order kernel diagrams denoted as H and F represent either a direct process 
(F) in which a particle-hole pair enters, interacts and exits the event, or an exchange 
process (H) in which a particle-hole pair enters, annihilates, and creates another pair 
which then exists the event. The second order kernel diagrams -81-3 and can be 

deciphered similarly. Diagram Bi describes a partially screened interaction between a 
particle and a hole, while diagrams B2 — B3 reflect the exchange symmetry, or describe 
processes in which the original particle-hole pair is annihilated, and another emerges 
the event. Diagrams Xi_3 relate only to partial exchange processes in which either an 
entering particle (^1,2) or hole (-B1.3) gets annihilated, and another particle or hole 
created during the process exits. 

We will now discuss the effect of the temporal structure of an approximation 
to the excitation spectrum which is given by the pole structure of the response 
function. We know that a time-local kernel only reproduce poles whose locations 
are given by differences in the renormalized, single-particle excited and ground state 
energies [HI [531 [131 [H]- The excitations corresponding to these poles are then 
called one-particle excitations, while all remaining ones we refer to as many-particle 
excitations. However, we can also view these one- to many-particle excitations in 
terms of some single-particle basis which has an appropriate energy label. In this case, 
excitations can have one- to many-particle character depending on their representation 
in this basis. In a wave function based theory, excitations are characterized with help 
of the basis expansion of the many-particle excited, and ground states. On the other 
hand, in many-body perturbation theory, we examine the diagrammatic structure of 
a response function in terms of Green's functions represented in this single-particle 
basis. 

In figurelHwe show some low-order diagrams which are obtained by expanding 2nd 
Born Green's functions in terms of Hartree-Fock Green's functions. These diagrams 
represent time-nonlocal processes which contain multiple, simultaneous particle-hole 
pairs, or in other words describe many-particle excitations in the Hartree-Fock 




(a) self-energy insertions (b) kernel contributions 

Figure 4: Some diagrams containing multiple simultaneous particle-hole pairs, whose 
origin lies either in self-energy insertions or kernel diagrams. A simple line represents 
a Hartree-Fock Green's function in terms of which the top and bottom rows describe 
two- and three-particle excitations, respectively. The direction of time is from left to 
right. 
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basis. Such additional excitations are clearly not present within the Hartree-Fock 
approximation and can, in general, arise either from self-energy insertions or directly 
from the kernel as shown in figures Hal and l4bl respectively. Moreover, self-consistency 
ensures that some classes of these diagrams representing many-particle excitations are 
taken into account to all orders of perturbation theory. If some of these diagrams are 
dominant in the perturbation expansion, we can expect many-particle excitations to 
appear in the excitation spectra. 

3. Model and test systems 

We consider A^-electron lattice systems which are based on the Pariser-Parr-Pople 
(PPP) model pil and are described with a Hamiltonian given by 

H{t) = h^J{t)clc,0 +UY, n^tn^l + "^^^(^^ ~ ^^^^^ " ^) ' (1^) 

where c^'^ annihilates (creates) an electron of spin a at site i, while hi = c\^c.ia and 
hi = Xla'^*" ^-""^ ^^^^ density operators. The one-body Hamiltonian is given by 

hij{t) -5(,,j)+5,jVi{t) (18) 

where the hopping elements <^{ij) are equal to 1 for the nearest-neighbor sites on the 
lattice and are zero otherwise. The symbol Vi{t) = v6i^iS{t) denotes a time-dependent 
potential. The two-body, electron-electron interaction is described with local and 
nonlocal interaction matrix elements denoted with U, known as Hubbard-JJ, and Wij, 
such that Wii = 0, respectively. These interaction matrix elements and the positive 
background potential Z are chosen to address two kinds of systems itemized below. 

(i) The Hubbard model is specified by choosing the off-diagonal interaction matrix 
elements, as well as the positive background potential, as zero {wij =0, Z = Q). 
The model represents a system with a short-range, or more precisely a contact 
two-body interaction. 

(ii) The PPP model, on the other hand, is introduced to address the effects of long- 
range interactions, which are described by the Coulomb type matrix elements 
Wij = U/2dij, where is a parameter describing the distance between sites i and j. 
We set the distance between nearest-neighbors on a lattice, as well as the positive 
background potential to one (c?(ij) = 1, Z = 1). 

We have, in particular, prepared lattice systems of two main types: linear chains 
and uniform rings, and conducted calculations by varying the interaction strengths, 
and number of electrons, but restricted ourselves to spin-compensated systems. In the 
next section, we summarize some of our results by considering a six site, six electron 
(1/2-fillcd) chain with long-range interactions, as well as a six site, two electron (1/6- 
fiUcd) ring with short-range interactions. Some essential features of both test systems 
are represented in figure [5l 

4. Results: addition, removal and excitation spectra 

We have chosen two finite lattice systems, described in section [3l to test the 
performance of the many-body approximations obtained from the Hartree-Fock and 
second Born self-energies. The observables used as a measure of the quality of 
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d chain 

u=i.o u=2.o ring 



2.58 - 
1.99 - 



3.36 
2.72 



[7=1.0 U = 2.Q 



a -0.11 ^ ^ „ ,1 

• • 

-1.58 

4 4i--1.36 -1.83 -fi- -1-67 

(a) spatial profile (b) Hartree-Fock energies 

Figure 5: The red sites in the spatial profiles of our test systems, shown in (a), denote 
the target site for a perturbation, while d = d^ij^ = 1 is the distance between nearest- 
neighbors. In (b), we show the Hartree-Fock single-particle energies for two interaction 
strengths {U = 1.0, 2.0). Notice that the energy spectra of the ring coincide with the 
non-interacting spectra up to a constant because of the discrete rotational symmetry. 



these approximations are the ground state energy, and the non-neutral and neutral 
excitation spectra calculated using the spectral and retarded response functions, 
respectively. The quality of our many-body schemes is assessed by comparing all 
approximate quantities against exact results. 

4.1. Numerical aspects 

All exact quantities are obtained by exact diagonalization (ED), or in other words 
calculated using the many-particle eigenstates of the system which are obtained by 
solving the time- (in) dependent Schrodinger equation. Our numerical implementation 
of exact diagonalization in these finite lattice systems relies on the well-known Lanczos 
routine ^5, 26, 27^. Approximate quantities, on the other hand, are obtained by 
calculating the equilibrium Green's function [161 HI] which yields the ground state 
energy, and from which the non-equilibrium Green's function is obtained by by time- 
propagation. These calculations are performed with a parallelized version 29 of the 
algorithm described in [301 131] ■ 

Some additional remarks on calculation of spectral properties are in order. These 
properties are, in practice, computed either by direct evaluation of addition, removal, 
and excitation energies (only in ED), or by finite-length time-propagation of a many- 
particle ground state (ED) and an equilibrium Green's function (MBPT). The latter 
approach leads to time-domain spectral and response functions which are transformed 
into their frequency-domain counterparts. As a side-effect of the finite propagation 
length, oscillatory noise-like features, as well as broadening of peaks, appear into these 
frequency-domain spectra. Such defects can be reduced by increasing the propagation 
length denoted with T and reached with a time step A. In figures, result related to the 
direct evaluation scheme will be denoted with ED-D, while time-propagated results 
are denoted with ED, HF or 2B. 
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4.2. Addition and removal energies 

The spectral function relates to so-called quasi-particle properties, in particular, the 
frequency-domain structure contains information about particle addition and removal 
energies, and the diagonal elements represent the local density of states. The Lehmann 
representation of the ground state spectral function is given by 



N+l\ 



/V-1 /V-1" 



5{u 



(19) 



where g^L' ^ (K 
strengths, and ^^^^ = 



c- 1$^+!^ 



rpN 



and = {^o\cl^\'^'^~'-) are the oscillator 

the addition and removal energies. Symbols \^^) 



and denote the kth eigenstate and -energy of a A^-particle system, respectively. 

The spectral functions of the 1/2-filled chain with different interaction strengths 
(U = 1.0,2.0) are shown in figure [B] Both many-body approximations capture the 
lowest addition and removal energies well, irrespective of the interaction strength. 
Moreover, they open up the band-gap adequately, and reproduce correctly the particle- 
hole symmetry, which although being in general broken by the long-range interaction, 
has been restored by the choice Z = 1 for the positive background. Differences between 
approximate solutions, as well between approximate and exact spectra, appear at 
higher energies, see addition energies A2-A5 or removal energies R2-R5, The 2nd 
Born approximation reproduces some of this additional structure present in the exact 
spectral function, while the Hartree-Fock method, merely capturing the single-particle 
energy spectra (Koopmans' theorem), fails to do so. 



3, 
<" 



R5i R4 

U = i1 .0 



R5 

u 



ED 
HF 



2.0 



A .;.ft....llt»^^iil 



Jj 



1 

CO 



I . 2B 



Figure 6: The spectral functions {Aii{uj) = Ai^^i^{uj) = Aii^ii{uj)) of the 

six site/electron (1/2-filled) chain with long-range (PPP) interactions. Exact 

addition/removal energies (ED-D) with a non-zero oscillator strength are given by 
the dashed, vertical lines. (T = 125, A = 0.02 - 0.025) 



The addition and removal energy spectra of the 1 / 6-filled, six site (2 electrons) 
ring are given in figure [T] Again, both approximations capture the lowest addition 
and removal energies well, but now only when U = 1.0, and reproduce correctly the 
broken particle-hole symmetry. As the interaction strength is increased, the Hartree- 
Fock approximation both underestimates the removal energies and overestimates the 
addition energies which, altogether, signals a failure to reproduce both the 2-electron 
ground state, and the lowest 3-electron states correctly. The ground state energies 
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given in table [Tb] verify these speculations, as well as the fact that both of our 
approximations describe the 1-particle system correctly, and therefore the removal 
energies describe solely the quality of the ground state. Although both approximations 
reproduce the addition energies (A1,A3) with the highest intensities, our correlated 
approximation gives more accurate results. However, only the 2nd Born approximation 
is able to quantitatively account for the additional, low-intensity structure which 
consists of a pole between high-intensity ones (A2) and two higher lying poles (A4,A5). 








-1 



■^^6!M^t^■ 

A3i, 



ill 



ED 
HF 
2B 



Figure 7: The spectral functions (Aii(w) = Ai-f^i^{uj) = Aii^n{uj)) of the six 
site, 1/6-filled (two electrons) ring with short-range (Hubbard) interactions. Exact 
addition/removal energies (ED-D) with a non-zero oscillator strength are given by the 
dashed, vertical lines. (T = 125, A = 0.02 - 0.025) 



Overall, we have observed that the Hartree-Fock approximation gives 
quantitatively incorrect high-energy spectra, which is simply due to fact that there 
are fewer single-particle energies than eigenenergies of the true interacting system. 
The second Born approximation, on the other hand, capable of accounting for the so- 
called correlation induced quasi-particles, does reproduce some of warranted additional 
structure, and additionally captures correctly shifting intensities as interaction 
strength is increased, see poles A2 /R2 of figure |6] or A2 of figure [71 In the next 
section, we will address the question of how well these fundamental differences affect 
the excitation structure evaluated at a higher S'E/SG-level. 



u 


So 


HF 


2B 


u 




HF 


2B 


1.0 
2.0 


-6.217 
5.540 


0.7% 
3.2% 


0.1% 
1.0% 


1.0 
2.0 


-3.867 
-3.782 


0.9% 
3.1% 


-0.1% 
-0.9% 



(a) 1/2-fillcd, six site, PPP chain (b) 1/6-fillcd, six site, Hubbard ring 



Table 1: The exact ground state energies Eq and relative errors Ehp/2B / Eq — 1^ where 
Ehf/2B is the ground state energy in the Hartree-Fock/2nd Born approximation. 
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4.3. Excitation energy spectra 

The neutral excitation energies can be accessed through the density(-density) response 
function, which is given by 

Xia.jfi (t, t') = Xfaia,jf}]l3 > (20) 

where the retarded response function is defined in equation In the fohowing, 

we consider time-translationahy invariant systems in which the response function 
depends only on a relative time, that is Xia,ji3it,t') = Xia,ji3{t — t'). The connection 
to the excitation energies becomes explicit in the Lehmann representation, in which 
the imaginary part of the ground state density response function can be written as 



Q'tn(xiaj73(t^)) = -TI"^ /n.ia/nj^ ^(5 (w - ri„) - 5 (cj + Vln)^ , 



(21) 



where fn,ia = ('^'o' I"-*" I'^n') is ^^'^ oscillator strength, and r2„ = — Eq the neutral 
excitation energies. Such equation is valid as long as the eigenstates of the A^-particle 
system, and hence the oscillator strengths, can be chosen real valued, which is the 
case for our time-reversal invariant systems. 

However, in practice, we calculate the ground state density response function 
Xijit) = J2ai3 Xia,jp(t) by mcans of time-propagation. First, we perturb the system 
with a spatially and temporally local potential Vi(t) = Vj6ijS{t), where j denotes the 
target site for the perturbation, v is the strength of the perturbation and 5(t) a Dirac 
delta function. Then, we use equation (1121) which states that the density response 
function is given by 

X^,{t)^':^^^^l^ + Oiv,), (22) 

where denotes the ground state site density, and n,;(t) = J^a^ia.iai^) ^i^^ 
density of the perturbed system evaluated at time t. Although nonlinear effects 
are always present, they are easily reduced by ensuring that the magnitude of the 
perturbation lies well in the linear response region. In practice this is verified by 
doubling the perturbation, and checking that also the response doubles to a sufficient 
accuracy. 

Before presenting neutral excitation spectra, we will introduce some auxiliary 
quantities which are useful for the discussion of the results. 

(i) The ground state energy can be used to judge whether a possible deviation from 
an exact excitation energy results from a failure to describe the ground, or excited 
state correctly. Although we cannot access the ground state energy at the / SG- 
level, a first order approximation denoted with Eqm is provided by the Galitski- 
Migdal functional ([8]) evaluated with the 2nd Born self-energy. 

(ii) We prove in [Appendix A| that any diagrammatically well-defined, self-consistent 
many-body approximation satisfies the f-sum rule, that is the density response 
function satisfies equations (I1.12p and (|1.13p . In our case, this simply means that 

'^ij — '^ij — '^ij J (23) 

where clj^*^"-* = J2ai3 '^iaj^^ ■ '-'^ ^^^'^ hand side of this equation, we have a 
commutator (see equation (|1.14p ) defined as 



c, 



^ ^jiTij + liihvj - Sij y2 [hjklk] + Ijkhkj) , (24) 



k 
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where hij and 7^- = J^ap 7*" j73 ^''^ ground state, one-body Hamiltonian and 
reduced density matrix, respectively. Whereas, on the right hand side, we have 
the first time-derivative and frequency- moment of the density response function 
which are given, respectively, by 

clf = -^tX^3it)\o+ ' (25 a) 

00 

c]f = J duj a;3m(xij (w)) . (256) 

We have also checked numerically, by fixing a method (ED, HF, 2B) and 
comparing the commutator ([M)) against the time-derivative (|25ap . that the f- 
sum rule is satisfied up to a numerical accuracy better than 0.01% for all of our 
methods (ED, HF, 2B). Moreover, since our response functions are obtained by 
finite length time-propagations, we have reserved the comparison between the 
commutator ([24| and frequency moment (j256p as a test of the completeness of 
our excitation spectra. 

As the frequency moment should equal to the commutator, we can use the latter 
to estimate how well an approximate method is able to reproduce the true spectra. 
In table [51 we compare exact (ED) and approximate (HF, 2B) ground state 
commutators (|24p . Such comparison shows that the second Born captures, in 
all cases, more precisely the value of the exact commutator, and is therefore 
expected to reproduce better quality spectra than the Hartree-Fock. 



u 


ci 


HF 


2B 


U 


ci 


HF 


2B 


1.0 
2.0 


1.755 
1.733 


1.1% 
3.9% 


0.4% 
2.1% 


1.0 
2.0 


1.324 
1.307 


0.7% 
2.1% 


-0.2% 
1.2% 



(a) 1/2-filled, six site, PPP chain (b) 1/6-filIed, six site, Hubbard ring 



Table 2: The relative error c^/c^ — 1, where cj^x = cl^y-^^ and X =HF or 2B. 



(iii) Lastly, we need some means to characterize excitations which appear in the 
exact excitation spectra. As the expectation value of an occupation number 
operator contains information about single-particle state occupations in a 
compact statistical format, it is a suitable tool for the task. Our occupation 
number operator is given by 

where d\^J annihilates (creates) an electron from (to) the single-particle eigenstate 
{i, a} of the Hartree-Fock Hamiltonian. We can obtain information about single- 
particle transitions by comparing the expectation values of this operator for 
the ground and excited states of the many-particle system. However, since 
excitation spectra can comprise degenerate excited states which have varying 
spectral intensities, we need a construct which can distinguish the character of an 
excitation in such a situation. For this purpose, we define the weighted occupation 
numbers 

ArHF(£;)^Tr(p,,(£;)nr), (27) 
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where Tr is the trace over a complete set of many-particle states, and 

Pt][^) = 7 7 ) (^^) 

L^k,Ek=E Jk,iajk,jp 

a properly normalized density operator which projects any state to a degenerate 
subspace, with eigenvalue i?, of the many-body eigenstates. The oscillator 
strength fk^ia is defined below equation (pT|). In practice, we calculate the 
differences 

ANf^liE) ^ Nf^liE) ~ Nf^^iEo) , (29) 

where Eq is the true ground state energy. This quantity tells us what kind of 
single-particle transitions are involved in the many-particle transitions from the 
ground state to the excited states with energy E. For example, consider a spin- 
compensated two-particle and -level system with nearly one-determinental ground 
state. Then weighted occupation numbers for the lowest single-particle state 
whose values are close to one or two describe one- or two-particle excitations, 
respectively. 

The neutral excitation spectrum of the 1/2-filled chain is shown in figure|S]for two 
{U — 1.0, 2.0) interaction strengths. The low interaction spectra (see the upper panel) 
are very similar, that is both many-body approximations reproduce the intensities and 
frequency structure of the exact response function quite accurately, although the 2nd 
Born approximation slightly better. As the interaction strength is increased (see the 
lower panel) both many-body approximations still reproduce quantitatively correct 
results, although some low-energy structure, see excitations El and E2, has started to 
shift towards the lower-end of the spectrum. This shift is a possible sign of a failure 
to describe the true ground state, as backed-up by the fact that the squared overlap 
of the Hartree-Fock ground state with the true ground state is only 0.86. However, 
since the correction Eqm — Eq ^ 0.06 is small, such shift cannot be attributed solely 




Figure 8: The imaginary part of the density response function xiii'-^) of ^ six 
site/electron (1/2-filled) chain with long-range (PPP) interactions, for which the f- 
sum rule is satisfied up to 99.9 — 100.2% when c\'i is used. Exact excitation energies 
with a non-zero oscillator strength (ED-D) are given by the dashed, vertical lines, 
(ui = 0.01, T = 150, A = 0.02 - 0.025) 
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as a breakdown of ground state. Overall, dominant features, like the intensities, and 
energies of low-energy excitations, as well as finer details, such as the diminishing 
intensity of the third dominant pole (E3) of the response function as a function 
of the interaction, are replicated more accurately by the correlated approximation, 
while the mean- field picture even fails to capture the subtleties. The fact that even 
the Hartree-Fock method reproduces all five, dominant excitation energies (E1-E5) 
correctly implies that we are dealing with one-particle excitations, as verified by the 
weighted occupation numbers of figure |9al We will shortly consider the question how 
do the many-body approximations cope with a more complicated excitation spectrum, 
in which additional excitations of many-particle nature arise. 



2 

m 1 



5 -1 
-2 



2 

UJ 1 


5 -1 
-2 



1 2 3 4 5 6 


1 2 3 4 5 6 


1 2 3 4 5 6 


1 2 3 4 5 6 


1 2 3 4 5 6 


E1 


E2 


E3 


E4 


E5 


(a) Excitations of the 1/2-fillod chain at [7 = 1.0, see figure[8] 
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1 ^1 — ^1 — ^1 ^1 1 


1 2 3 4 5 6 


E1 


J E2 


□-^ 

E3 


U E4 


J E5 



(b) Excitations of the 1/6-filled ring at U = 1.0, see figure [TOl 

Figure 9: The relative, weighted occupation numbers ANk{E) = AN^^i^{E) are shown 
for each Hartree-Fock state {k = 1, —6), and excitation E1-E5. The histograms whose 
filled (positive) or depleted (negative) parts sum up to 0.5 — 1.5 or 1.5 — 2.0 describe 
a state with a singly- or doubly-excited character, respectively. 



We give an example of such a system by introducing a less rigid, but energetically 
more simple, 1/6-filled, six site ring, whose neutral excitation spectrum is given 
in figure [10] for two {U = 1.0,2.0) interaction strengths. The results portray, 
irrespective of the interaction, an adequate agreement between the exact and second 
Born excitation spectra, but on the contrary, an utter failure of the Hartree-Fock 
approximation to account correctly for the structure of the response function. Two 
excitations, one at a lower energy (E2) and another at a higher (E4-E5), which are 
completely absent in Hartree-Fock, indicate the presence of two-particle excitations, 
which cannot be described with a time-local approximation. The weighted occupation 
numbers, shown in figure [Qb] confirm these speculations: only two excitations. El and 
E3, relate to singly-excited states, while all others, in particular E2,E4 and E5, have 
a strong doubly-excited state character. Therefore, the highest excitation of Hartree- 
Fock attempts to mimic a practically non-existent lowest to highest state transition, 
see figure I5bl while the 2nd Born approximation successfully reproduces the two- 
particle excitations, the lowest one being quite accurate, while the higher ones are 
shifted towards lower energies. As a disadvantage of the correlated scheme, a possibly 
spurious excitation appears on left-hand side of E3 when U — 2.0. Such non-physical 
defects have also been seen in other recent studies [131 E] , where they are associated 
with a response kernel which does not have a proper mathematical structure. 
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E 4 



U = 1.0 



U = 2.0 



ED 
HF 

2B 



1 2 3 4 5 6 

(0 

Figure 10: The imaginary part of the density response function xii(a;) of a six site, 
1/6-filled (two electron) ring with short-range (Hubbard) interactions, for which the 
f-sum rule is satisfied up to 99.8 — 100.1% when c}f is used. Exact excitation energies 
with a non-zero oscillator strength (ED-D) are given by the dashed, vertical lines. 
{vi ^ 0.01; T = 150, A = 0.02 - 0.025) 



5. Summary and conclusions 

We have studied the performance of the time-dependent Hartree-Fock and second 
Born approximations by investigating the spectral properties of some finite, strongly 
correlated, lattice systems. Our purpose has been two-fold: firstly, to test 
approximations used in transient quantum transport, and secondly to gain insight into 
excitations in MBPT, both using the Kadanoff-Baym equations. We have calculated 
ground state energies, as well as spectral and density response functions with our 
approximate methods, and compared these against exact results, which were obtained 
by exact diagonalization. 

The results show that both approximations perform well in a simple half-filled 
system. That is, they reproduce correctly the true excitation spectra, although the 
quasi-particle properties are not replicated as accurately, especially in the Hartree- 
Fock approximation. However, at a lower filling, we observe two-particle excitations, 
which by construction cannot be captured with the Hartree-Fock method, but are 
generated correctly in the second Born approximation. Such a difference is also seen in 
the spectral function in which additional quasi-particles are formed in this correlated 
approximation. Overall, the second Born approximation is consistently in a better 
agreement with the exact results, although signs of spurious excitations, also in more 
asymmetric and strongly interacting systems, call for further research. 

We conclude that a time-local approximation can be inadequate for highly 
correlated lattice systems at a low-filling, which should have an impact on quantum 
transport, as particle number on device varies. On the contrary, already the 
simplest time-nonlocal, self-consistent approximation performs more reliably. More 
complicated, but computationally realizable, many-body approximations, as well as 
investigation of the effects of self-consistency, remain as motives for a future exposition. 
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Appendix A. Thomas-Reiche-Kuhn/f-sum rule 

The exact density(-density) response function satisfies a consistency relation known 
in atomic physics as Thomas-Reiche-Kuhn (TRK) sum rule, or in condensed matter 
physics as the f-sum rule. It is known [32] that the TRK/f-sum rule holds in the 
spin-position basis for self-consistent MBPT. However, there does not exist a proof in 
the literature that it is satisfied in a lattice system, in which locality gets redefined 
and the notion of electromagnetic gauge invariance becomes less transparent. In the 
following, we prove that a response function obtained using equation ()12|) fulfills the 
TRK/f-sum rule also in a lattice system. 

The electromagnetic gauge freedom is described in real-space by the 
transformation {v{r,t),A{r,t)) {v{r,t) — dtA{r,t), A{r,t) + VA(r,t)), where 
{v{r,t), A{r,t)) and A{r,t) are the gauge potential and function, respectively. We 
can extend this transformation to a lattice system by introducing a new gauge 
function Aj(z), where j is a collective site and spin index. This function is defined 
on the Keldysh contour of figure [T] on which it satisfies the boundary conditions 
Ai(to) = Ai{to — i/?). Now, we define the gauge transformation h{z) — > h-^{z) where 

h'^iz) = ^e'(^'(-)-'^^(-))/i,,(z)clc, -^9,A,(z)ri, (A.l) 

ij i 

is the gauge transformed one-body part of the Hamiltonian. Moreover, using the 
equation of motion (|4]), we can show that 

G^-(z,z') = e^'^'''>G,,{z,z')e-'^^^''^ . (A.2) 

Here, we found essential that the interaction matrix elements are two-index quantities, 
since then the exponential phase factors cxp(iAi(z)) cancel at each interaction vertex. 
This cancellation guarantees that 

(z, z') = e'^'(") (S[G]) (z, z')e-'^^("') , (A.3) 

which in practice is only true if the one-body Green's function is calculated self- 
consistently. The Hamiltonian (jA.l[) is given, to first order with respect to Ai(z), 

by 

h^{z)^h{z) + v{z) + --- (A.4) 
where we defined the one-body potential 

Wy(z) = i/ijj(z)(A,(z) - Aj{z)) - Sijd^A^iz) . (A.5) 
Consequently, the first order variation of the one-body Green's function can be written 
as 

SGfj{z,z') = ^ Jdz Lij^ik{zz',z)vki{z) 

7 



kl 



= 1 

kl 



7 



dz \^Li.i^ki{zz',z)hik{z) - Lijjk{zz',z)hki{z) 
iSkidzLijjk{zz', z)^ Ai(z) , (A.6) 
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where Lij ki{zz', z), is the generahzed response function defined in equation (I13p . The 
second hue above was obtained by using partial integration, as well as the boundary 
conditions Ly^/cj (zz', io — i/?) = Lij^u{zz' .t^) and Ai(to — i/3) = Ai(to)- We can also 
expand equation (jA.2[) to first order with respect to Ai(z), and equate the result 
against the variation (jA.6[) . This leads to a result given by 

{6,k5{z' ,z") - 8jk5{z,z"))G,,{z,z') 

= ^ (jiki{z")Lijjk{zz', z") - Lij^kiizz', z")hik{z")j - idz" Lij^izz' , z") , (A.7) 
I 

which is nothing but the generalized Ward identity relating the generalized response 
function to the one-body Green's function. At the limit z' = z^ using equation (jl4p . 
the Ward identity can be split into three real-time equations, one for each Keldysh 
component. These equations are given by 

^dt'xi,,{t,t') = ('^'^•Ki')x|a.(i,^') - xl,,{t,t')hi,{t')) , (1.8&) 

where ^ij{t) and xfj kii^^^') denote the one-body reduced density matrix and 
greater/lesser component of the real-time response function which are defined in 
equation ^ and p.ip . respectively. Consequently, the density response function 

x^,M = e{t - t'){xl,,{t,t') - xl,j{tX)) (1.9) 

satisfies an equation given by 

i^t'X^Jit,t') = 5{t t'){xl^^{t,t') - xl,,{t,t')) + 9{t t') ^ {h,,{t'){xl,^{t',t) 

k 

-Xt.,k,(t'^t)) ixhkit'^t) - Xl,kit'^t))hk,{t')) , (1.10) 

where we used the Ward identity (|1.86p . Moreover, evaluating this equation at equal 
times t' = t — rj with the limit 77 0+ taken afterwards and using the symmetry 
Xl,ki^t,t') = xtl,^J{t',t) (see equation ^) leads to 

i5t'Xy(^,i')L,^t- 

= E {hok{t){xl,uit,t) ~ xfkM^,t)) - (xr,-,,(t,i) - X%At^t))hk,{t)) . (1.11) 

k 

If the initially unperturbed system is time-translationally invariant then the density 
response function depends only on a relative time, or in other words Xij (^: t') = 
Xij{t ^ t')- Then the Ward identity (|1.8ap allows us to write the Thomas- Reiche- 
Kuhn sum, or the f-sum, rule in time-domain as 

- <9tX»i(0|t=o+ " '^i^^^J + ^i'^^'i ~ E i^jklkj + Ijkhkj) ■ (1.12) 

k 

Furthermore, by assuming that the density response function is analytic in the upper 
half of the complex plane, the TRK/f-sum rule can be written in the frequency- 
domain [32] as 

j du! a;3m(xij(w)) = /ijiTy + 7ji^y - ^ij E i^Jklkj + Ijkhkj) ■ (1-13) 
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Note that in the text, we refer to the right hand side of these equations as the ground 
state commutator cj^ since in an exact theory it can be written [3 3) as 

<), (1.14) 

where \^o) is the iV-particle ground state, hi the site density operator, and H{t) the 
many-particle Hamiltonian. 
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